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1 Introduction 

In this paper we demonstrate how a self-consistent treatment of the physics of the early 
universe (from temperature T ~ 10 MeV down to T ~ 0.1 eV) enables cosmic microwave 
background (CMB) observations, in concert with light-element abundances, to be used as 
new probes of beyond-standard-model (BSM) physics in the neutrino sector. In a sense, 
these probes are tantamount to probes of the Ci^B (relic Cosmic Neutrino Background). 
Recent observations [1-4] of the CMB and observationally-inferred primordial abundances of 
deuterium and helium [5-8] formed in big bang nucleosynthesis (BBN) already place tight 
constraints on both the cosmological standard model (CSM) and BSM physics. However, 
future observations will usher in a higher level of precision with even better prospects for 
BSM probes, see for example Ref. [9] . 

We anticipate that observations will bring about an overdetermined situation where 
BSM physics may manifest itself if it is present. The existing or anticipated observations 
and measurements of greatest utility for the present purposes are: (1) high-precision mea¬ 
surements of the baryon-to-photon ratio, or the equivalent baryon density {ujh = Hf,/i^); 
(2) high-precision measurements of the “effective number of relativistic degrees of freedom” 
(Neff); (3) high-precision measurements of the primordial deuterium abundance (D/H) from 
quasar absorption lines; (4) measurements of the primordial helium abundance (Lp) directly 
from CMB polarization data; and (5) measurements of the sum of the light neutrino masses 
(Y^niu), i.e. the collisionless damping scale associated with neutrinos. 
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The physics that determines the relic neutrino energy spectra in weak decoupling and 
that of primordial nucleosynthesis affects observables of the CMB. These distinct, disparate 
epochs, however, depend not only on the values of parameters describing the cosmology (such 
as ojb, iVgff, relevant cosmic constituents, etc.) but also on the parameters derived from these 
base, input parameters (such as primordial abundances, recombination history, etc.). This 
is the basis for what we will term “self-consistency.” The requirement for self-consistency 
between the BBN and CMB epochs, for example, in current data analyses depends on the 
parametrization of the energy density in terms of A^eff upon the baryon-to-photon ratio. 
These two parameters, as measured at photon decoupling, are the sole determinants of the 
primordial helium abundance at the epoch of alpha particle formation (T ~ 0.1 MeV) in 
“standard” (i.e. zero lepton numbers, no BSM physics) BBN calculations. Though this 
procedure is adequate for the CSM and standard model physics, here we argue that it can 
be insufficient to probe varieties of BSM physics when confronting model cosmologies with 
next-generation, high-precision CMB and light element abundance data. In the case of the 
helium and self-consistency mentioned above, the relationship between the helium yield 
in BBN and Ngg is a non-trivial function of the interplay of expansion rate and neutron-to- 
proton ratio, as is well known [10]. The latter ratio is, in turn, a sensitive function of the 
zze and zzg energy distribution functions, and these can be affected by BSM issues like lepton 
numbers, flavor mixing, sterile neutrino states, heavy particle decay, etc. 

Handling and analyzing the observed data while imposing self-consistency over multiple 
epochs in the early universe can require new procedures. The approach developed in Ref. 
[10] and employed in Planck XVI [11] has been highly successful. There, self-consistency is 
imposed approximately by the addition of a term in the log-likelihood function. A very small 
theoretical uncertainty ensures the posterior distributions are close to the corresponding 
theoretical values. However, when we impose self-consistency among the different epochs 
the need may arise to solve for the various derived parameters iteratively. Returning to the 
helium/Agff example discussed above, Yp and the recombination history of the universe both 
depend on the radiation energy density, usually parametrized by Aeff- Subsequently, the 
recombination history is affected by the the primordial abundance of helium. The values 
of these two derived parameters, in turn, affect CMB observables, which recommends an 
iterative and therefore self-consistent approach. 

Ideally, we would want a procedure to self-consistently treat neutrino and BSM physics 
from the post-QCD epoch to the onset of non-linearities in large scale structure (LSS). Re¬ 
cently, the quantum-kinetic equations (QKEs) governing neutrino flavor evolution in dense 
environments have been derived from first principles [12-14]. Such a program to treat neu¬ 
trino physics in the early universe would incorporate a QKE treatment through weak decou¬ 
pling at the very least. Weak freeze-out and BBN necessitate a neutrino-energy Boltzmann- 
equation or QKE treatment fully coupled to a nuclear reaction network. The recombination, 
photon decoupling, and advent of LSS epochs require a Boltzmann-equation treatment of 
neutrino clustering. Verification of any models related to neutrinos and BSM physics would 
then rely on agreement with direct cosmological observables, including but not limited to: 
primordial abundances; CMB power spectra; and the total matter power spectrum. 

As outlined, this is a challenging undertaking. We therefore employ a limited approach 
to show why self-consistency is required and how such a treatment can be efficacious in prob¬ 
ing some issues in neutrino sector physics. To that end, we calculate the ratio of sound horizon 
to photon diffusion length, r^/r^, as a simple parametrization of the CMB, as described in 
detail in Sec. 3. It would, of course, be preferable to compute the full CMB power spectra in 
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the self-consistent manner described above. There are two limiting factors that temper this 
ambitious proposal, however. The first is that the observable rs/r^ is largely insensitive to 
the poorly constrained equation of state of the dark, vacuum energy component while the 
CMB power spectra are not. Additionally, current computations of the CMB power spectra 
[15] would need to be generalized to the BSM scenarios we contemplate here. 

Our procedure for calculating derived cosmological quantities utilizes the neutrino oc¬ 
cupation probabilities. Since we are developing the capacity to compute the effects of BSM 
physics that couples to the active neutrino sector, we do not restrict the form of the neutrino 
distribution functions to that of equilibrium Fermi-Dirac distributions. General forms are 
permitted and may be handled analytically or numerically. Neutrino occupation probabilities 
are taken in the present work in the flavor eigenbasis during weak decoupling, weak freeze-out 
and BBN.^ During recombination and last scattering, we transform the occupation proba¬ 
bilities to the mass eigenbasis. When we consider the case of massive neutrinos, the statistic 
for the sum of the light neutrino masses, ^ my, is simply the sum of the three lightest mass 
eigenstates. 

We might further clarify the fact that this limited approach does not rely on the as¬ 
sumption that the radiation energy density can be described in terms of the single parameter 
Aeff. Indeed, an original motivation for using r^/r^ as a proxy for Aefj (see Sec. 3) is the fact 
that the usual definition for Aefj does not apply to non-equilibrium neutrino distributions. 
A corollary of this observation is the possibility that may depend on the scale factor 
a{t)] that is, we do not assume that Agfr is independent of scale factor. When extended, our 
approach is able to handle BSM scenarios such as the decay of massive sterile neutrinos and 
other weakly interacting massive particles, which result in neutrino distributions that may 
differ substantially from a Fermi-Dirac distribution [16]. 

Our approach aims to be general; it does not rely on particular cosmological models but 
is appropriate to the class of A cold-dark matter (ACDM) models. It extends and explicates 
recent work [9, 10, 17, 18] on the consistent incorporation of precision observations of the 
CMB [11] and observations of the light element abundances of helium and deuterium (mass 
fraction Yp and relative abundance D/H). The deuterium abundance [6, 8] is now more 
precisely measured by a significant factor (4 or 5) in relative precision than Yp. Our results 
indicate that deuterium is sufficiently well determined to be incorporated into CSM analyses 
as a fixed prior. 

The present work is a first result in an ongoing campaign to incorporate neutrino energy 
transport into a self-consistent treatment of BBN and recombination. Our approach is being 
implemented in Fortran90/95 as a suite of codes under the working title of “BBN Unitary 
Recombination Self-consistent Transport” (burst), which will be made publicly available for 
use on parallel computing platforms using openmpi. The current work is a prerequisite in 
an ongoing collaboration to develop codes that consistently handle BBN and neutrino energy 
transport from weak decoupling to the advent of LSS. 

The paper is outlined as follows. In the following section. Sec. 2, we discuss our treat¬ 
ment of BBN and relevant CMB physics. Section 3 explains the self-consistent, iterative 
approach to determining Aeg using r^/r^. Section 4 employs a dark-radiation model to test 
the accuracy of predictions made by burst against observations. Section 5 investigates a 
number of ‘test problems’ from BSM physics suited for burst. We conclude in Sec. 6 with 
future applications of the results given in this work. Throughout this paper, we use natural 

^This is a matter that is properly resolved by utilizing a QKE approach [12], the subject of a future study. 
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units where h = c = ks = 1- 


2 Background 


2.1 Big-bang nucleosynthesis 

Recently there has been a dramatic increase in precision in the determination of D/H [8] . This 
improvement in observational precision drives the need to improve the standard tools used 
to calculate observables during the BBN and CMB epochs. This development allows us to 
consider improved descriptions of nuclear physics and non-standard particle and cosmological 
models. 

We have developed a generalized BBN subroutine, part of the larger burst code. The 
current BBN routine handles general distribution functions for particles out of equilibrium. It 
numerically integrates the binned neutrino occupation probabilities to obtain thermodynamic 
variables and number and energy densities for Ue, Ug) W- R is worthwhile, therefore, 

to summarize our approach to BBN even though it is substantially similar in some respects 
to that given in Refs. [19-21]. 

We assume that the cosmic fluid is homogeneous and isotropic and that the comoving 
number density of baryons is covariantly conserved. Neutrinos can experience, in addi¬ 
tion to gravity, essentially arbitrary interactions within the Boltzmann transport approx¬ 
imation.^ We allow for the possibility that neutrinos may decouple from the plasma with 
non-equilibrium distributions. This assumption implies that there may be deviations from the 
Fermi-Dirac momentum distribution [22-24]. In fact, it was this observation that provided 
initial motivation for developing the current approach. 

BBN codes evolve light nuclide abundances Yi{t), defined as 




ni{t) 

nb{ty 


( 2 . 1 ) 


as a function of comoving time t in the background of the Friedmann-Lemaitre-Robertson- 
Walker (FLRW) geometry. Here ni{t) is the proper number density of nuclide i = n, p, ^H, 
^He, ^He,... and nfe(t) is the proper baryon number density. We will focus on two specific 
nuclides in this paper: the '^He mass fraction {Yp = and the relative abundance 

(D/H = Yp/Yff). The evolution starts from a time when the plasma temperature T is near 
30 MeV. Weak equilibrium obtains at this temperature. The time dependence of the metric 
is determined by the energy density p{a) as 


= I = \/§^' 

where H(a) is the Hubble parameter, rupi is the Planck mass, and the ‘dot’ indicates differ¬ 
entiation with respect to the FLRW-coordinate time t. The energy density is computed as 
the integral of the single-particle energy over the momentum distribution: 

p{a) = J a)^p^ + mj (2.3) 

^However, the code can handle limited generalizations of the Boltzmann approximation to incorporate 
effects associated with neutrino quantum kinetics. 
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where the sum over j reflects all species contributing to the energy density, including but 
not limited to: photons, baryons, dark matter, e^, and neutrinos. Consistent with the 
above discussion, we do not need to assume a well defined temperature for any of the cosmic 
species. The time dependence of the distribution function /i(p; a) is indicated by the presence 
of the scale factor a{t) in its argument. Evolution in time of the distribution functions is 
accomplished by solving transport equations, such as the Boltzmann equation, in FLRW 
geometry. 

Given initial conditions for the temperature T and momentum distributions of the cos¬ 
mological constituents (assumed to be in equilibrium through weak and strong interactions), 
the BBN code determines the evolution, with respect to the scale factor a{t) (related to time 
as dt = da/{aH{a))), of the temperature of the plasma T{a), the electron chemical potential 
/ie(a), and the nuclide abundances relative to baryon number Yi{a). 

The interactions of the light nuclear species is governed by the nuclear reaction network.^ 
The reaction network, which is determined by a chosen set of nuclides and the thermally 
averaged reactivities na^na 2 {(^jiaVa), is proportional to the rate of change of the number 
density of nuclides participating either in the initial state n^- or final state where i,j 
indexes particles (up to 3) in the initial a or final /? channel for the process a ^ P; that is, 
a and /3 are two- or three-body reaction channels and is the nuclide of the channel a. 
The nuclear reaction network includes nuclides with mass number ^4 < 9. Our code allows 
for the inclusion of additional nuclides with ^ > 9, but we maintain the smaller network 
as the larger network provides no new insights for this paper. The nuclides are taken to 
be in thermal equilibrium with the photon-electron plasma. We are currently employing an 
updated [25] version of the reaction network from Ref. [21]. We also couple in all relevant 
and neutrino-induced weak interactions (charged and neutral current) [26]. 

2.2 Cosmic Microwave Background 

The physics of BBN and photon decoupling are related through the epoch of recombination 
by several mechanisms including, importantly. He i and He ii recombination. The recombi¬ 
nation rates depend sensitively on the amount of ^He created in BBN [10, 27, 28]. The ^He 
dependence also determines the photon diffusion damping wave number or, equivalently, 
the photon diffusion length, defined as = vr/Zcrf. This quantity depends strongly on the 
recombination history through the free-electron fraction, Xf,{a), which, in turn, depends on 
the helium abundance (see Sec. 2.2.3 below). 

We relate the diffusion length directly to an observed angular size 6d through the angular 
diameter distance Da at the epoch of last scattering (as described in detail below). Since 
Da depends on the vacuum (dark) energy equation of state, which is poorly understood, we 
employ the ratio OslOd, equivalent to Vg/rd (see Eqs. (2.6) and (2.13)), where Og is the angle 
subtended by the sound horizon, which eliminates explicit dependence on the dark energy 
equation of state. 

Modern computations of the recombination history of the early universe are available in 
accurate and well tested codes such as HyRec [29], CosmoRec [30], and the fast, approxi¬ 
mate computations given in Ref. [31] by the code RecFast [32]. Our recombination history 
compares well with RecFast, agreeing at the level of < 3% over the recombination epoch. 
The present accuracy suffices for the purposes of the current exploratory work. 

^The nuclear reaction network should be derived from reaction cross sections that are governed by the 
principles of quantum mechanics, such as unitarity. We are currently developing a unitary reaction network 
for application in future work. 
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2.2.1 Sound horizon 

The sound horizon is defined as 


r^'yd 

iio,'yd) — / dcL ■ 
Jo 


a‘^H{a)^3{l + R{a)) 
where a^d is the scale factor at photon decoupling."^ The ratio R{a) is given as 


(2.4) 


R{a) 


3pb(a) 

4.p-y{a)' 


(2.5) 


where pi, and p^ are the baryon rest mass and photon field energy densities, respectively. 

The angular size of the sound horizon , determined from the spacing of the acoustic 
peaks in the CMB temperature power spectrum, is related to the sound horizon by the 
angular diameter distance DAia^yd) at photon decoupling (scale factor a^d) as 


Bg {CLryd) 


Cl'yd 


rsia-fd) 
Da (o!"yd) 


since the angle Bg is small. The angular diameter distance is 


T'a(o) 



da' [a'^H{a')]-\ 


( 2 . 6 ) 


(2.7) 


The quantity DA{a^d) depends on the vacuum (dark) energy equation of state, which is not 
very well understood. Our approach will eliminate dependence on this poorly constrained 
component of the energy density. 


2.2.2 Free-electron fraction 

We have written an independent code to calculate the free-electron fraction, W. The results 
agree well with Ref. [32] (RecFast). In fact, the agreement is within 2% for most of the range 
^ ^ 10“^ (10"^ ^ z ^ 10^). The code allows significant deviations from the model 

parameters of ACDM; it should be used with caution, however, since the effective three-level 
treatments for helium and hydrogen recombination have been optimized for near-standard 
model parameters. 

The number density of free and total electrons is 

=np + UHe „ + 2nHe (2.8) 

(2.9) 

where nb, Up, nueu, ^^Hem are the proper number densities for baryons, protons, singly- and 
doubly-ionized helium, respectively. We write the free-electron fraction as 

(free) 

Tt' ' 

^ -(tot)" ^ II + 2WHe III, (2T0) 

Tl^ 

so defined to take values in the range 0 < W < 1. 

^Reference [11] takes a^d —t a* where the optical depth is unity. 

^Reference [11] takes 9s —t 6^*. 
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We follow Refs. [33, 34] and consider the simplification of the multi-level hydrogen and 
helium atoms to that of an effective three-level system which includes the ground n = 1 
state, the first excited n = 2 states, and the continuum. All other excited states are assumed 
to be in equilibrium with the 2s state. We treat He ii recombination approximately [31] 
(via the Saha equation) since it is essentially complete at the advent of the epoch of He i 
recombination. The He lll contribution, therefore, in the Boltzmann equation for He ll is 
negligible. 

Boltzmann equations for H ii and He ii contain a thermally-averaged cross section 
and relative velocity (cru). We use Case B recombination coefficients for the {av) of both 
neutral hydrogen [35] and helium [36]. Along with a Saha equation for He iii, the Boltzmann 
equations for H ii and He ii are a coupled set of ordinary differential equations constituting 
a recombination network to model the ionization history of the universe prior to photon 
decoupling. 


2.2.3 Photon diffusion damping length 

In the tightly coupled limit, photon diffusion damping is characterized through the damping 
wave number [37-40] 

2 ^ da 1 R^{a) + ^{1 + R{a)) 

Jq a?‘H{a) ane{a)(TT 6(1-|-i?(a))^ ’ 

where ctt is the Thomson cross section. Here, we have assumed that moments of the temper¬ 
ature fluctuation higher than the quadrupole make a negligible contribution in the linearized 
Boltzmann equation for the photon distribution. 

Equation (2.11) requires the free-electron fraction, ne(a), determined in the recombina¬ 
tion history, as discussed in the previous section. The free-electron fraction, over the course of 
the recombination history, depends strongly on the primordial helium mass fraction Yp. The 
fraction Yp, in turn, depends on the relativistic energy density, which is often parametrized 
in terms of the parameter 


Pr — 2 





( 2 . 12 ) 


Equation (2.12) is only valid after the epoch of annihilation. The point that BBN and 
recombination are related is well known [10] but has not been implemented self-consistently 
as a constraint for general, BSM physics model cosmologies. We return to it after describing 
the relation between and kd to directly observable quantities given by the CMB power 
spectrum. 

The observed diffusion angle Od{a^d) [H] is related to the diffusion damping length, 
Td = tt /kd as 

0d(aT,d) =aT,rf-^^, (2.13) 

with the same stipulation regarding the smallness of the angle in Eq.(2.6). 


3 Tg/vd as a proxy for A^eff 

In this section we describe our method for determining A^efr in detail. We introduce two 
variants of A^eff- When referring to the radiation energy density equation (2.12), which 
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takes as an input the quantity A^eff; we designate A^eff as Aigg*^\ When considering general 
cosmologies, perhaps with BSM physics, we deduce the value of Aigff from the observable 
quantity = Osl^d-, described in this section, and designate it as iVeff. The simplest 

cosmologies for which Eq. (2.12) obtains, having negligible neutrino mass, standard model 
constituents and no energy transfer between species have A^gg*^^ =fVeff. 

We consider a test input cosmology that is non-standard yet substantively similar to 
ACDM. We proceed by determining Yp at temperature T ~ 0.1 MeV using the BBN network 
of BURST [25]. The principal observational cosmological input at this time is ujh - Also input 
and incorporated into the BBN subroutine is any neutrino and BSM physics that constitute 
the test cosmology. Subsequently, we compute the recombination history of the universe from 
early times (a/ao ~ 10“^) to the current epoch. Specific observational inputs include w;,, cug 
(where for cold dark matter ojc = Yp, and Hq (the Hubble constant, Hq = H{ao)). 

For the purposes of the present discussion, other inputs of particular importance include 
neutrino occupation probabilities (as output from a neutrino transport calculation of weak 
decoupling that is fully coupled to BBN) and neutrino rest masses. The neutrino energy 
density of the neutrino seas is calculated by writing the occupation probabilities in the mass 
eigenbasis [41]. We emphasize the fact that is not input as a base parameter; this is 

of paramount import in the present approach. The recombination history, ne(o) determines 
the optical depth as a function of scale factor: 

rao 

'^(o) = / —p^ a'ne{a')aT (3.1) 

Ja a 

We define the scale factor at photon decoupling such that T(ajd) = 1- In this definition, 
we do not include the effects of cosmic reionization when calculating ne(o) for use in Eq. 
(3.1) [11]. We apply a^d and the input cosmology to equations (2.4) and (2.11) to compute 
the sound horizon and photon diffusion length, respectively. We arrive in this way at a ratio 
for our input cosmology. 

Our immediate objective is to determine a value of Ngff (here termed Neflf) corresponding 
to this value for (rs/rrf)^“P^. We map out a range of values of r^/r^ that correspond to the 
same input cosmology as that used in calculating with one significant difference. 

We parametrize all of the neutrino and BSM physics into the single A^gg^^ parameter. We 
then use A^gg^^ to calculate the radiation energy density in Eq. (2.12) to determine the Hubble 
rate. We vary A^gg^^ to compute the function 

rs/rd = rs/rdluJb,u}c,Yp,.. •; Ng^g*"^], (3.2) 

shown in Fig. 1. Since r^/r^ is a one-to-one function of we may invert Eq. (3.2) to 

obtain The final step is to evaluate the previous function with our 

input cosmology ratio, i.e. Ngg = N^j^^rg/vd = (Ts/rrf)(“P)], to obtain a value of N^s- As 
an example, we take the best-fit values from Ref. [11] combined with WMAP Polarization 
data (1000^ = 1.04136 & 1000^ = 0.161375) to obtain (r^/rrf)('°P) = lOO^^/lOO^rf = 6.45304. 
This corresponds to a value A^eg = 3.31 on Fig. 1, in line with the best-fit value Aeg = 3.25 
(3.511^:?^ at 95% limits) [11]. We again note that for the simplest cosmologies the two 

generally distinct functions [r^/r^] and Nes['>’s/i"d\ reduce to the same function and have 

'vi?' = N.S. 
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Figure 1. Ratio of the comoving coordinate of the sound horizon radius rg to that of the photon 
diffusion length as a function of for cosmological parameter values Yp = 0.2425, Wf, = 0.22068, 
and LOc = 0.12029. 


Figure 1, which shows the function demonstrates an important constraint 

between phenomena occurring during the epochs of BBN and recombination/photon decou¬ 
pling. For a given input cosmology {i0b,ujci ■ ■ ■), the graph of r^/rd as a function of 
requires a computation of rig (a) to obtain rd for each value of 

We can understand Fig. 1 qualitatively by a simple scaling argument. We expect that, 
as the radiation energy density increases with increasing the sound horizon and the 

diffusion length will decrease with the increasing Hubble rate H[a). The sound horizon 
decreases due to the increased energy density driving a more rapid expansion and a decrease 
in the sound speed, Cg = [3(1 + The diffusion length increases, naively, due to 

a decrease in the scattering rate driven by the reduced Hubble time. Caution should be 
taken when using such naive scaling arguments. For example, the non-trivial dependence of 
the recombination history leads to counterintuitive effects in the Wff dependence on ^ rrii, 
[42]. Ref. [42] demonstrates that where a naive scaling argument would suggest an increase 
in Wff with increasing the non-trivial dependence of the recombination history on 

^ mi, implies decreases monotonically and rapidly with increasing ^ (See Sec. 5.1). 

4 Verification of cosmological parameters with dark radiation 

We adopt the point of view advocated in Refs. [10, 27, 28] that the constraint provided by 
predictions of BBN should be incorporated simultaneously with constraints due to recombi¬ 
nation effects in the extraction of cosmological parameters. We also require, as previously 
discussed, that BBN and recombination be solved iteratively. In this section, we demonstrate 
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the self-consistent extraction of these parameters by employing a simple model for the radi¬ 
ation energy density that avoids solving, for example, the Boltzmann equation, for the set of 
distribution functions of the cosmic constituents, in particular the neutrinos. 

The model we explore in this section is identical to ACDM except for one additional 
constituent: dark radiation. The dark radiation energy density is radiation at all epochs 
and does not interact with the other energy-density constituents through any force except 
gravitation. 

The total energy density as given in terms of radiation, matter, and vacuum energy 
components is 


p{a) = pr{a) + pm{a) + pv{a), (4.1) 

which depend on the scale factor a as a“^, a~^, a°, respectively, as long as there is no energy 
transfer between the species. We assume for the purposes of this section that the neutrinos 
are massless, always acting as radiation energy density. The matter energy density consists 
of contributions from baryons and cold dark matter. Modeling the matter as a pressureless 
gas and observing that the comoving matter energy density is conserved, we write the proper 
energy density as 


3H^ / an \ 3 

Pm= ° + . (4.2) 

The vacuum energy is the least understood of the energy densities. Assuming the 
universe to be critically closed, the sum of the three energy densities must be equal to 
the critical energy density, specified only by the Hubble rate at the current epoch. Hence, 
Pv = Pc — Pr,o — Pm,o- The vacuum energy density is negligible at all epochs of interest in this 
paper but is included for completeness. 

The radiation component is given as: 


Pr — P'y + Pv + Pdr 


(4.3) 


where p-y is the photon energy density, p^ is the neutrino energy density, and pdr is the 
dark-radiation energy density. We parametrize pdr as 




30 


(4.4) 


where (idr is the dark-radiation parameter and is always assumed to be non-negative. In 
principle, we could entertain negative values of 5dr since it is an adjustable parameter of pr- 
Such a change requires a fundamental reworking of the ACDM model so that 
These non-standard cosmologies obtain when considering, for example, neutrino oscillations. 
These models, however, are not continuously connected with our model at ddr = 0, for any 
values of the parameters, thus motivating the maintenance of ddr > 0. 

We write = 3 -|- and assume that the contribution to from pi, is 3. 

Then the contribution from pdr is given as We see then that AA^^^^ = ddr- This is 

simply a restatement of the fact from the previous section, Sec. 3, that NeS = 3 -|- AA^g^^ = 
Agg^^ for ‘standard’ cosmologies. It is, therefore, unnecessary for this simple dark-radiation 
model, to deduce NeS from r^/r^ since A^g^^ =Aeff by construction. This model of dark 
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radiation is the usual model applied, for example in Ref. [11] and we explore, in this section, 
the predictions of the present burst code to verify our results against those of prior results 
within the community. We will use A^eff; for the remainder of this section, to denote the 
“effective number of relativistic degrees of freedom®.” 

Figures 2-5 show the results of computations in which the four parameters ujb, Yp, D/H, 
and iVeff are varied. We begin by varying the two model inputs: uJb and Ngg. The upper 
panel in Fig. 2 shows the dependence of A^gff on uJb for curves of constant Yp; the vertical 
band is the Planck value of cob = 0.02207 ±0.00033. The figure is generated by first choosing 
a value for the baryon density in the range 0.004 < u)b < 0.029. Each selected value of the 
dark radiation parameter in the range 3 < A^eff ± 4.5 allows for the prediction of Yp and D/H 
by parametrizing the radiation energy density as in Eq. (2.12). The values so obtained are 
plotted as contours in the N^s-^^b plane in the upper and lower panels of Eig. 2. The solid curve 
is the preferred value Yp = 0.2465 ± 0.0097 of Ref. [7], which is a selection of observations of 
metal poor extragalactic H ii regions. The contours are spaced by roughly 0.0097/3 showing 
that A^eff is not strongly constrained by values of Yp alone; this is a manifestation of the 
degeneracy of A^eff and Yp. For example, at Yp = 0.2465 ± 0.0035, corresponding to the 
contours closest to Yp = 0.2465, the range allowed Algff is nearly consistent with both the 
standard, calculated value A^eff = 3.046 and the Ref. [11] derived value of A^eff = 3.30 ± 0.27. 

Predictions of the primordial deuterium abundance are a much more sensitive constraint 
upon allowed values of Ngg. This can be seen in the lower panel of Fig. 2. The solid line 
contour with 10® x D/H = 2.530 ± 0.04 corresponds to the recent measurement of Ref. [8]. 
Contours in this hgure are separated by the one standard deviation of Ref. [8] . There are two 
points of interest regarding the deuterium figure. First, as noted in Ref. [8], observation of the 
primordial component of deuterium is precise enough to begin to constrain the microscopic 
physics of the thermally averaged nuclear reaction rates and their cross sections. Additionally, 
given the precision of the current and forthcoming deuterium measurements and the strong 
dependence of on its value (at constant uJb), we advocate using D/H as a prior, over Yp, 
for future base model parameter searches. 

The degeneracy between N^g and Yp is again evident in Fig. 3. Each plot explores the 
D/H vs. LOb contour space, where the upper plot contains contours of constant Yp and the 
lower plot contains contours of constant N^g. The shaded bands in each figure indicate the 
one-sigma observations of ojb and D/H from Refs. [8, 11], respectively. Deuterium is not an 
input parameter into our model. We compute it by choosing a baryon number and iteratively 
change the dark-radiation parameter, (5dr until matching the chosen deuterium target. The 
outputs from the process are Nf,g and Yp. Values of oJb, D/H and Neg are in satisfactory 
agreement with the standard cosmology at the precision of current observations. 

The quantities D/H and LOb are the tightest observationally-constrained parameters we 
are currently investigating. Eigure 4 shows two plots in the Yp-N^g plane with contours of 
constant cob (upper plot), and contours of constant 10® x D/H (lower plot). The horizontal 
band in each figure indicates the one-sigma observation of Yp from Ref. [7]. Like D/H, Yp is 
not an input into our model. Consequently, we adopt the identical iterative method for Yp in 
Eig. 4 as we do for D/H in Fig. 3. For the oJb (upper) plot of Fig. 4, the solid contour line is 
the best-fit value of Ref. [11] with nine-sigma spacing of the contours. The contours exist in a 
subspace of the Yp-N^g plane which is well within current observations, but nevertheless could 
span a range of radically different physics. Similarly, the bottom plot shows the 10® x D/H 

®We refer to N^g as the effective number of relativistic degrees of freedom although there are factors that 
complicate this interpretation, among them the temperature parameter, and fermionic nature of neutrinos. 
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Figure 2. (Top) iVgg plotted against Wf, for contours of constant values of Yp (labeled by mass 
fraction). The solid curve is the preferred value of Ref. [7]. The contours are spaced by AYp « 0.003. 
(Bottom) Ngg versus Wf, for contours of constant values of 10^ x D/H. The solid curve is the preferred 
value of Ref. [8]. The contours are spaced by A(10^ x D/H) = 0.04. In each case, abundances are 
determined in a self-consistent BBN calculation. 
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UJb 


Figure 3. (Top) 10® x D/H plotted against uJb for contours of constant Yp. The solid curve is the 
preferred value of Ref. [7]. The contours are spaced by AYp « 0.003. (Bottom) 10® x D/H versus Wf, 
for contours of constant iVgff. The contours are spaced by AiVeff = 0.3 


value of Ref. [8] as the solid contour with the other contours spaced hfteen-sigma apart. 
Clearly, Yp and Agfr do not constrain the cosmological model as tightly as ujf, and D/H. This 
observation indicates the import of using the next generation of 30-meter class telescopes and 
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Figure 4. (Top) Yp plotted against iVgff for contours of constant uJb- The solid curve is the best-fit 
value of Ref. [11]. The contours are spaced by Aojt = 0.003. (Bottom) Yp versus Ngg for contours of 
constant 10® x D/H. The solid curve is the preferred value of Ref. [8]. The contours are spaced by 
A(10® X D/H) = 0.6. 


CMB observation to better determine the light element abundances, particularly, Yp with 
high precision. 
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Ub 


Figure 5. Yp plotted against w;, for contours of constant iYgff. The contours are spaced by 
ATVeff = 0.33. 


Figure 5 shows how changes in the Yp-ujb plane. The shaded bands in each figure 
indicate the one-sigma observations of Ub and Yp from Refs. [7, 11], respectively. We may 
conclude from this plot that there could exist many different values of Ngg consistent with 
the observations of Yp and cuf,. However, we caution against such a conclusion without 
considering the effect on D/H. In fact, we choose not to include a figure of contours of D/H 
in the Yp-ojb plane because the strong sensitivity of D/H to ujb produces too large a range of 
values for Yp to be useful as a constraint of the cosmological model. 

All calculations show a consistency between w;,, AeSj D/H, and Yp to a conservative limit 
of two-sigma error range in each observation. We expect the uncertainties in each observation 
to improve in the coming years with large ground-based CMB experiments [43, 44] and 30- 
meter class telescopes [45-47]. Future high-precision measurements may result in tensions 
for the best-fit values of Ub, iVeg, and D/H. These tensions could be indications of the need 
for more precise theoretical and numerical approaches or could signal the presence of physics 
beyond the standard model. As it stands here, the bottom panels of Figs. 2 and 3 shows 
that the tension between uJb and D/H cannot be resolved with the addition of extra radiation 
energy density. Uncertainties in nuclear reactions may produce disagreement between Ub 
and D/H, allowing D/H to become a probe of nuclear physics. It is also possible that 
the spectroscopic determination of D/H may be subject to small systematic errors only 
recognizable at such precision. An exciting prospect is the need to revise the CSM to resolve 
tensions with the observations of D/H and cufe, possibly leading to the conclusion of BSM 
physics active during BBN. 
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Figure 6. Determination of plotted against Wb at constant AiVeg. The contours are spaced by 

« 0.01 in values of AiVgg. All contours correspond to AiVgg < 0. 

5 Examples for neutrino sector BSM physics 

We describe three examples of unresolved issues in BSM/CSM physics, which require the fully 
self-consistent parameter determination described in previous sections. We consider, in turn, 
models incorporating neutrino rest mass, sterile neutrinos, and non-zero lepton numbers. We 
use, throughout this section, the observationally-inferred definition of Aigg, Agff- 

5.1 Neutrino rest mass 

Section 4 details a self-consistent treatment of the BBN observables Yp, D/H, and 
The sum of the light neutrino masses, denoted Y^rriu, has no bearing on the determination 
of primordial abundances in BBN calculations due to the high temperatures relevant there. 
Here, however, we explore the epochs and energy scales in the history of the universe as¬ 
sociated with the rrii, energy scale in order to investigate the relationship between Y 
and the other four observables of interest (cuf,, iVgg, Yp and D/H). Specific examples of such 
epochs that we might consider include the surface of last scattering (z ~ 1100) and the 
advent of LSS (z < 10). We focus on the surface of last scattering and implications for the 
CMB in this paper. 

Reference [42] (hereafter GFKPI) investigates the effect of neutrino rest mass on Ngg 
using the burst suite of codes. Conventional estimates based on the energy density added by 
non-zero neutrino rest masses suggest an increase in However, using the method out¬ 

lined in Sec. 3, GFKPI shows a decrease in Ngg. The decrease is due to an effect on the recom¬ 
bination history stemming from an increase in the Hubble rate, which results in a larger free- 
electron fraction. This counterintuitive result is termed the “neutrino-mass/recombination 
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(z^MR) effect.” The z/MR effect manifests itseff only in a self-consistent treatment, such as 
that employed by GFKPI. 

In addition, GFKPI investigates the dependence of the z^MR effect on ujh- We revisit 
this physics here in preparation for a discussion on the effect of non-zero lepton number Li, 
later, in Sec. 5.3.2. Increasing cjf, leads to an enhancement of the z/MR effect when ^ rrii, 
is held constant, as is evident by the curvature of the contours in Fig. 6. The enhancement 
is a consequence of the effect that changing ujh has on the the recombination history. We 
might naively expect a larger change in the Hubble rate relative to the massless neutrino 
case resulting in an enhanced z/MR effect for the smaller ujb case. This is opposite to that 
observed in Fig. 6. This result also is counterintuitive based on expectations from a simple 
scaling of the energy density and the resulting change in the recombination history [42]. 

The origin of the enhancement of the z^MR effect can be understood by considering 
a simplification of the Boltzmann equation that determines the recombination history [Eq. 
(2.10)]. We take Yp = 0 for the purposes of this argument since the z/MR enhancement 
is insensitive to Yp, as we have verified numerically for the ranges of parameters we are 
considering. In this simple scenario, Xg = Xp and we obtain an expression for the change in 
the free-electron fraction: 

^ = (l-Xe)/3-X2n(*°t)«(2), (5.1) 

where f5 = a^^'^irrieT is the ionization coefficient and is the recombination 
coefficient with AQ = 13.6 eV. 

Equation (5.1) for the recombination history shows that the free-electron disappearance 
rate is proportional to the total electron number density, which in turn is proportional to 
ujb through Eq. (2.9). Equation (2.9) also shows how relates to Yp. Note that Ub and 

Yp affect differently. However, due to the relative insensitivity of Yp to uib, the Ub 

dependence dominates in Eq. (2.9). The increase in energy density from ^ 0 and the 

increase in the free-electron disappearance rate combine to alter the recombination history 
so as to enhance the z/MR effect for increasing cob- 

5.2 Sterile neutrinos 

We next consider the possibility that there exists either single or multiple sterile-neutrino 
species, which could have profound implications in cosmology. We entertain two possibilities 
of either light or heavy sterile neutrinos. 

5.2.1 Light sterile neutrinos 

Observations of neutrino events in large scintillating detectors may have revealed anomalies 
that could be interpreted as sterile neutrinos with rest masses ~ 1 eV [48-50]. We 
investigate the presence of a single sterile neutrino in the early universe by employing a model 
where the sterile state populates a thermal Eermi-Dirac shaped distribution with temperature 
parameter Tg, possibly through flavor mixing. The sterile neutrino temperature, Tg is taken 
to be less than or equal to the active neutrino temperature T^. The ratio Tg/T^ is assumed 
to be the same throughout weak decoupling, BBN, and recombination. For this analysis, 
we do not investigate smaller active-sterile neutrino mixing angles with resultant non Fermi- 
Dirac-shaped energy spectra [41, 51]. Future work will consider such physics [52]. 

Figure 7 displays contours of constant Neg. The vertical axis is the ratio Tg/T^, and the 
horizontal axis the sterile neutrino rest mass rriy^. We maintain the ratio Ty/T = (4/11)^/^, 
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Figure 7. Ratio of the sterile to active neutrino temperatures, Ts/T^, plotted against for 
contours of constant iVeff for ^ rriy = 0.06 eV. Horizontal dotted lines show the prediction if the 
sterile neutrino was massless, i.e. = 0 


assuming covariant conservation of entropy, starting at the end of the epoch of e^-annihilation 
and continuing throughout the remainder of the history of the universe. The dotted lines 
show the expectation from the dark radiation analysis of Sec. 4 without employing the self- 
consistent, iterative approach developed in Sec. 3. The deviation of the contours from the 
dotted lines is again due to an effect similar to the zzMR effect but, in this instance, due 
to the sterile state. Fig. 7 takes the sum of the active neutrino masses to be 0.06 eV. This 
is inconsequential for large Tg/Ty < 1. For Tg/Ty < 0.1, AiVeg < 0 due to the zzMR effect 
in the active neutrino sector. As a consequence, the contour for = 3 is not coincident 
with the axis. Since rriy^ is too small to be of any significant kinematic effect during 
BBN, we need only compute BBN once for a given value of Tg/T^. During recombination, 
is kinematically important and affects the Hubble rate. Hence, for every point in the 
Tg/Ty-niy^ plane of Fig. 7 we calculate recombination. This figure clearly emphasizes the 
need for a self-consistent treatment between BBN and recombination when considering this 
BSM physics. 

5.2.2 Heavy sterile neutrinos 

Heavy sterile neutrinos that decay out of equilibrium in the early universe can affect weak 
decoupling and, as a consequence, primordial nucleosynthesis [16, 53, 54]. Sterile neutrinos 
in the rest mass range 0.1 GeV < < 1.0 GeV, with lifetimes ^ 1 s decaying during the 

weak decoupling, weak freeze-out, and/or BBN epochs can have constrainable, sometimes 
dramatic, cosmological effects. 
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Such sterile neutrinos have mass and vacuum mixings with constrained by 

accelerator and other laboratory oscillation experiments/observations [55-62], beta-decay 
experiments [63], and cosmological considerations, including constraints on and A^efr 

[64-69]. In fact, stringent constraints can be obtained from A^eff limits alone [16], as sterile 
neutrinos decaying out of equilibrium can lead to dilution (entropy production) which, in the 
weak decoupling epoch, can lead to distortions in the relic neutrino energy spectrum, affecting 
'^rriy, and have significant impact on the relativistic energy content and, hence, A 

sophisticated theoretical and computational treatment of dilution physics is a challenging 
endeavor. Even though these heavy sterile neutrinos may decay away before an epoch where 
T ~ 10 keV, they can nevertheless alter the relationship between Yp, D/H, and 

cuft, necessitating the need for a self-consistent treatment between the weak decoupling, weak 
freeze-out, BBN, recombination, photon decoupling, and advent of LSS epochs [70]. 

5.3 Lepton numbers 

We examine how lepton numbers affect the primordial abundances and iVefj. We define the 
lepton number Ly for a neutrino species p in a flavor eigenstate as 

_ Uu-np 

Ly — - 

n-y 

where Uy is the number density of neutrino species v, rip is the number density of anti¬ 
neutrino species u, and is the number density of photons. Here, for illustrative purposes, 
we take Ly^ = Ly^ = Ly^ = Ly. In fact, flavor oscillations will bring the various lepton 
numbers into near equality [71-73]. 

We use the comoving-invariant neutrino degeneracy parameter ^y = ^yjTy, where ^y is 
the chemical potential of neutrino v, to compute Ly. The present model assumes that the 
lepton number evolves through the epoch of annihilation only in response to the relative 
increase in n.y. We do not consider any BSM physics which could alter the difference riy — np; 
that is, we fix ^y throughout weak decoupling, BBN, and photon decoupling. We relate the 
degeneracy parameter to the lepton number using the following expression [51]: 

where C(3) ~ 1.202. The factor 4/11 in Eq.(5.3) implies that our lepton numbers refer to the 
post annihilation epoch, where Ty/T = (4/11)^/^. 


(5.3) 


(5.2) 


5.3.1 Effect on nucleosynthesis 

The helium mass fraction is sensitive to the neutron-to-proton ratio, n/p. We determine n/p 
by calculating the weak rates associated with neutrino-nucleon reactions, namely: 

Ue + n ^ e~ + p'^ 

Ue + ^ + n 

and in addition, neutron and inverse neutron decay: 

n ^ e~ + Ve + P~^ (5-6) 

The net rates in reactions (5.4) and (5.5) fall below the Hubble rate at the epoch of weak 
freeze-out. Weak freeze-out largely precedes the alpha-particle formation process in BBN, 


(5.4) 

(5.5) 
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Figure 8. Lepton asymmetry [Eq.(5.2)] plotted against ujb for contours of constant Ip. The 
contours are spaced by AYp = 0.015. 


though unlike the brief time/temperature range of a-formation, weak freeze-out occurs over 
several Hubble times at this epoch. The rates in reactions (5.4) through (5.6) are sensitive to 
the neutrino and distributions. We follow Ref. [21] to evolve T and the electron chemical 
potential in order to maintain equilibrium between the electrons, positrons and photons. For 
the electron-flavor neutrinos, we use the comoving invariants aT^y and to compute the 
neutrino distributions. We set = 0 as neutrinos of sub-eV rest mass remain ultra- 

relativistic throughout weak freeze-out. 

Figures 8 and 9 show the helium mass fraction and the relative deuterium abundance, 
respectively. Each plot is in the Ly-un, plane for contours of constant primordial abundance. 
The relationships between lepton number and nucleosynthesis are well known [74, 75]. In¬ 
creasing leads to an overabundance of neutrinos compared to anti-neutrinos. The forward 
rate of reaction (5.4) freezes-out after the forward rate of reaction (5.5). The imbalance low¬ 
ers n/p which lowers Yp as seen in Fig. 8. The decrease in n/p also leads to a decrease in 
D/H, although deuterium is not as sensitive to L^y as helium. However, D/H is known to 
much higher precision than is Yp. 

Comparing with recent observations [7, 8], the two light element abundances achieve 
consistency at 2(T. Yp prefers a value of Lj/ < 0 whereas D/H prefers a positive value of Ly. If 
future observations of the light-element abundances were to show a larger disagreement than 
2(7, lepton numbers of identical value could not solely rectify the tension. Future analyses will 
consider scenarios with multiple facets of BSM physics including non-zero lepton numbers 
[76]. This analysis will use A^eff as a discriminating factor. 
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Figure 9. Lepton asymmetry Ly plotted against Wf, for contours of constant 10® x D/H. The solid 
curve is the preferred value of Ref. [8]. The contours are spaced by A(10® x D/H) = 0.04. 


5.3.2 Effect on A^gff 


We consider how two aspects of non-CSM/BSM physics ^ 0 and/or ^ / 0) modify 

Ngg. If we set the nentrino rest mass to 7^ 0, we can investigate whether the z/MR 

effect still applies with non-zero L,y. Figure 10 shows the changes to Agfr in the '^rrii,- 
cob plane for three values of Li, = —0.05,0,0.05 corresponding to blue, magenta, and red 
contours, respectively. The non-zero lepton number increases Ngg for small values of ^ nii^ 
in accordance with Ref. [9]. However, for values of ^ nii, ~ 1.0 eV, the i^MR effect overwhelms 
the extra energy density from more particles to lower below three. 

Figure 10 also shows that, despite the total energy density being insensitive to the sign 
of L^, the contours for non-zero values of Li, with opposite sign do not overlap because Yp 
depends sensitively on its value. Yp is largest for the blue contours, so more helium suppresses 
the pMR effect. 

Note that taking Ly ^ Q conflates the interpretation that the effect of ^ rriy is identical 
to perturbations in the matter power spectrum, which are used in calculating the suppression 
of power on small scales. This is borne out by the present model where the ^ rriy statistic 
cannot be equated to the cosmological measurement. In our model, ^ rriy is simply the sum 
of the active vacuum neutrino mass eigenvalues. The observationally determined value of 
^ rriy depends on quantities other than the sum of the active neutrino masses such as their 
energy distributions. 
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Figure 10. plotted against Wf, for contours of constant AiVeff. The blue contours are for 

Li, = —0.05. The magenta contours are for L^, = 0. The red contours are for = 0.05. Solid 
contours are for positive values; dashed contours are for negative values. 

6 Conclusions and outlook 

Cosmological considerations are a key route to exploring BSM physics. This is especially true 
for the neutrino sector, where there are many outstanding questions and where laboratory 
experiments are limited in what aspects of this physics can be addressed. In this paper we 
have argued that a self consistent treatment of BSM issues, across all epochs from weak 
decoupling to photon decoupling, is the best way to take advantage of the expected coming 
increase in precision of CMB measurements and observationally-inferred primordial abun¬ 
dances of the light elements. We employ a limited prescription to link the salient features of 
self consistency between early-time neutrino dynamics and the surface of last scattering. We 
couple the weak decoupling and nucleosynthesis of early times to CMB observables, including 
baryon-to-photon ratio (equivalently Wb), sound horizon, and photon diffusion length. 

We have shown that such a self consistent treatment is necessary, in part because new 
neutrino physics can alter the relationships between different cosmological epochs. For ex¬ 
ample, cjb and other CMB observables affect the calculated yields of deuterium and helium. 
In addition, the calculated relic neutrino energy spectra after weak decoupling affects the 
predicted value of iVefj at photon decoupling. The principal tool in our analysis is the suite 
of BURST codes for nucleosynthesis and neutrino interactions and energy transport. 

A case in point is our investigation of the relationship between neutrino rest masses, i.e. 
and the four potential observables ujh, iVeff) Ip, and D/H. This analysis reveals the 
“neutrino-mass/recombination” (uMR) effect first described in Ref. [42]. The uMR effect is 
below the threshold of current CMB capabilities, but may not be in future observations [77] . 
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There have been spectacular advances in the measurements/observations of neutrino 
properties. We know the neutrino mass-squared differences and three of the four parameters 
in the unitary transformation between the energy eigenstates (mass states) and the weak 
interaction eigenstates (flavor states) of the three active neutrinos (only the CP-violating 
phase remains unmeasured). As active neutrinos mix in vacuum and have non-zero rest 
masses, the question indubitably arises of whether there exist “sterile” neutrino states. If 
indeed sterile neutrinos do exist, we acknowledge that the parameter space of mass, vacuum 
mixing angle, and number is enormous. However, sterile neutrinos could have profound 
effects in all of the epochs under study in this paper. This possibility makes a self consistent 
treatment of these effects a powerful basis for constraining sterile neutrino states. 

In this paper we have considered scenarios for both “light” (mass ~ 1 eV) and “heavy” 
(mass ~ 0.1 — 1 GeV) sterile neutrinos. In the former case we consider cases where the sterile 
neutrino relic energy spectra are Fermi-Dirac black-body shaped, though with a temperature 
parameter Tg differing from that characterizing the relic energy spectra of active neutrino 
species. We show here that the z/MR effect has interesting consequences and that this case 
demands a self consistent treatment of recombination and BBN. Additionally, heavy sterile 
neutrino decay out of equilibrium can lead to dilution and high energy relic active neu¬ 
trinos, and both of these features potentially can have dramatic and constrainable effects 
on CMB-epoch observables. This implies that CMB observations can indirectly probe the 
Cz/B and explore active-sterile mass/vacuum-mixing parameter space unavailable to current 
accelerator-based experiments. 

We have also studied the effects of non-zero lepton numbers on the relationship between 
CMB observables, nucleosynthesis, and neutrino physics. Our conclusion is that the primor¬ 
dial deuterium abundance is a potentially powerful probe of lepton number. However, an 
eventual CMB-only measurement of the primordial helium abundance Yp will be the most 
powerful probe of lepton number and many other issues. Determining Yp from CMB ob¬ 
servables will require a sophisticated self-consistent approach to BBN, neutrino physics, and 
photon decoupling transport physics. 

Finally, our study has revealed a potential tension between Aefr, oob, and the primordial 
deuterium abundance, D/H, inferred from high redshift QSO absorption systems. In fact, 
if the advent of 30-m class telescopes in the near future allows for a decrease in errors in 
observationally-inferred D/H to the ~ 1% level, while observed ojb and A"efr maintain their 
respective current central values, then tension is unavoidable. This may signal BSM or non- 
CSM physics, likely in the neutrino sector, or it could point to not understanding systematics 
in the damped Lyman-a cloud measurements of the isotope-shifted hydrogen absorption 
lines. We advocate using future instruments to explore the rich physics of weak decoupling, 
nucleosynthesis, and photon decoupling to discover what role BSM neutrino physics has in 
these epochs. 
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